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Figure 1: The first-order momentum field induced by a symmetric Gaussian 

bump. We took m = 1, A = 1, a = 1, P x = 2, P y = 0. (a) The 

G\ momentum field, (b) The G2 momentum field. (The 

perturbation is probably too strong for a first-order result to be accurate, 

but the purpose here is only to show what the first-order field looks like). 

Figure 2: Plot of H(6 X , 9 y , Jjjp-, Jjjp-) at various times for the 
Pullen-Edmonds Hamiltonian. (a) t = 0.0. (b) t = 0.1. 
(c) t = 0.2. (d) t = 0.3. 

Figure 3: Plot of H(6 x ,6 y , Jgn, Jgn) for the Pullen-Edmonds Hamiltonian 
for t = 00 in the first-order limit. 

Figure 4: Plot of H(9 X , 6 y , Jjjp, ^j-) at various times for the 
Pullen-Edmonds Hamiltonian with the (2, —2) resonance removed, 
(a) t = 0.0. (b) t = 0.1. (c) t = 0.2. (d) t = 0.3. 
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Abstract 

This paper presents a PDE-based approach to finding an optimal canonical basis 
with which to represent a nearly integrable Hamiltonian. The idea behind the method 
is to continuously deform the initial canonical basis in such a way that the depen- 
dence of the Hamiltonian on the canonical position of the final basis is minimized. 
The final basis incorporates as much of the classical dynamics as possible into an in- 
tegrable Hamiltonian, leaving a much smaller non-integrable component than in the 
initial representation. With this approach it is also possible to construct the semiclas- 
sical wavefunctions corresponding to the final canonical basis. This optimized basis 
is potentially useful in quantum calculations, both as a way to minimize the required 
size of basis sets, and as a way to provide physical insight by isolating those effects 
resulting from integrable dynamics. 
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1 Introduction 

Suppose we are given a nearly integrable Hamiltonian i?(q, p), where (q, p) represents some 
canonical representation of phase space (not necessarily ordinary position and momentum). 
We wish to find another canonical representation (Q, P) in which H is as close as possible 
to being integrable. Our motivation for this is twofold, and is connected to semiclassical 
quantum mechanics: First of all, from a purely numerical perspective, a representation in 
which H is as close as possible to being integrable leads to an optimized semiclassical basis 
with which to perform quantum calculations. The reason for this is that we can associate 
with P a quantum state |P). If we write H(Q, P) = H^(P) + # {1) (Q, P), then it is clear 
that is diagonal in the {|P)} basis, so that only is available to couple the various 
basis states. The coupling is given semiclassically by [g], 



which is as small as possible will minimize the couplings in the corresponding {|P)} 
basis, and thus will minimize the size of the basis required to perform a given calculation to 
some desired accuracy. 

Secondly, a canonical represenation in which H is as close as possible to being integrable 
provides physical insight. By incorporating as much of the classical dynamics as possible 
into an integrable Hamiltonian, this optimal canonical representation can help to isolate 
those classical and quantum effects resulting from integrable dynamics from those that do 
not. Thus this representation can isolate non-integrable quantum and classical effects such 
as dynamical tunneling and Arnol'd diffusion, respectively. Furthermore, we can also visu- 





where 




Fourier component at P t, F of A representation in 



2 



alize the quantum manifestation of the integrable classical dynamics via the semiclassical 
prescription for constructing a wavefunction: Given the generating function «9(q, P) from 
the initial basis (q, p) to the final basis (Q, P), we have, up to normalization, 

(q|P> Hdet^l^ exp[iS//l] (2) 

where det jf-§p is the well-known Van Vleck determinant. 

In the case of action-angle variables, the optimized representation of a nearly integrable 
Hamiltonian is termed an Intrinsic Resonance Representation (IRR), a term coined by Car- 
ioli, Heller, and Moller (CHM) [p], §]• In 1997 they published a paper detailing an algorithm 
for the construction of such a representation. The idea behind the CHM algorithm is to 
eliminate all the non-resonant terms of the Hamiltonian via an appropriate canonical trans- 
formation. This canonical transformation is obtained via a modified Chapman, Garrett, and 
Miller (CGM) method JIL 0, E|, which is essentially a Newton- Raphson scheme to find the 
invariant tori with a desired set of actions for a nearly integrable system. The remaining 
resonant and near-resonant terms are then re-expressed in the new basis. It is impossible 
to reduce the angle dependence any further, since this would result in the formation of 
resonance zones, which prevents a global action-angle description of the Hamiltonian. 

In 2001 the author published an alternative algorithm for finding the IRR basis [I]. The 
method introduced is a PDE-based approach which continuously deforms the initial action- 
angle basis in such a way that the angle dependence of the Hamiltonian is continuously 
reduced. It amounts to a gradient-descent algorithm in the limit of a first-order perturba- 
tion, and was therefore called the GDA method. Formally, the method does not distinguish 
between resonant, nearly resonant, and non-resonant terms, that is, the evolution is per- 
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formed on the entire Hamiltonian without any terms neglected. However, the evolution is 
such that the more non-resonant a term, the more strongly it is affected by the evolution. 
Thus, the non-resonant terms of the Hamiltonian are essentially killed off, the nearly reso- 
nant terms are reduced somewhat, while the resonant terms are essentially unaffected. The 
GDA method was successfully applied in Ref. f to 2, 4, and 6 degree-of-freedom systems. 

The GDA method circumvents two main drawbacks of the CHM method. First, the CHM 
method requires an a priori decision as to which terms are resonant and non-resonant. This 
leads to an ambiguity in the case of near-resonances. It could happen that a given term in 
the Hamiltonian must be considered resonant in order to get the Newton-Raphson scheme to 
converge. This leads to a somewhat artificial cutoff criterion, since a nearly resonant Fourier 
component should in principle still be reduced as much as possible, though not necessarily 
completely. Thus, unless the Hamiltonian has a few exact or near-resonances, it is not clear 
that the CHM method will give the optimized torus basis. 

Second, the CHM method requires the numerical evaluation of multidimensional integrals, 
and the numerical inversion of a nonlinear angle map at every iteration step. These numerical 
calculations slow the algorithm down. In contrast, the numerical calculations required by 
the GDA method are much simpler, so we believe that the GDA method is faster than the 
CHM approach (though in fairness it should be added that no direct speed comparisons have 
been made to date). 

Despite its advantages, the GDA method is limited to systems describable by action-angle 
variables. Furthermore, while the method generates the representation of the Hamiltonian 
in the final action-angle basis, it does not give the overall generating function S(q, P), trans- 
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forming from the initial (q, p) basis to the final (Q, P) basis. Thus, we could compute the 
energy spectrum arising from an optimized invariant torus basis, but not the corresponding 
semiclassical wavefunctions "living" on the IRR tori. 

This paper generalizes the GDA method to arbitrary canonical representations of phase 
space, and extends it to include the determination of the overall generating function S^q, P). 
It is organized as follows: In Section 2 we derive the generic evolution equation for the 
Hamiltonian starting from an arbitrary canonical basis (q, p). We then go on to derive 
the corresponding evolution equation for the overall generating function, with which it is 
possible to construct semiclassical wavefunctions. In Section 3 we consider the case of a 
nearly integrable Hamiltonian, and obtain the specific form our evolution is to take if we 
want to optimize the canonical basis used to represent the Hamiltonian. In Section 4 we 
consider the first-order limit of our PDE approach. In particular, we obtain a generalized 
first-order classical perturbation theory which coincides with classical perturbation theory 
in the case of action-angle variables. In Section 5 we consider some analytical examples. We 
look at free propagation in one dimension perturbed by a localized potential, and also in 
two dimensions perturbed by a symmetric Gaussian bump. We continue in Section 6 with a 
numerical example. Finally, we conclude in Section 7 with a summary of our results and a 
discussion of future research plans. 

2 The Evolution Equations 

In this section we shall derive the basic evolution equations for the Hamiltonian H and the 
overall generating function S. The idea is as follows: We start with an initial set of canonical 
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coordinates (q, p), which denote any global representation of phase space, and do not nec- 
essarily refer to ordinary position and momentum, respectively (non-global representations 
can also work, as long as we remain well within the region of phase space they describe. An 
example of this is action-angle variables for a system which can dissociate. Action-angle vari- 
ables should work as a valid representation as long as we remain well below any dissociation 
threshold). We continuously deform this system via a series of infinitesimal generating func- 
tions. The result is that our canonical representation is evolving with time, and is denoted 
by (Q t , at time t. As our canonical pair evolves, the functional dependence of H on the 
canonical pair changes. In addition, the overall generating function S(q, P; t) connecting the 
initial (q, p) to the current (Q t ,Pj) evolves as well. In subsection 2.1, we derive the PDE 
governing the evolution of H generated by this phase space deformation. In subsection 2.2, 
we derive the PDE governing the evolution of S. 

2.1 Evolution of the Hamiltonian 

Consider an arbitrary set of canonical coordinates. At time t, we're at system (Q t ,P t ). At 
time t+dt, we're at system (Q t+dt , P t+dt ). These are connected by an infinitesimal generating 
function F(Q t , P t+dt ) = Qt • Pt+dt + dtG(Q t , P t+dt ). Therefore, 



Qt+dt = Qt + dtX7pG(Q t , Pt+dt) 



(3) 



P t = Pt+dt + dtV Q G(Q t ,P t+dt ) 



(4) 



so to first-order, we obtain, 



Qt — Qt+dt — dtVpG (Qt+dt, Pt+dt) 



(5) 
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P t = Pt+dt + dtV Q G(Q t+dt , P t+dt ) (6) 
The Hamiltonian at time t + dt is therefore related to the Hamiltonian at time t via, 

H^ +dt \Q t+dt ,P t+dt ) = H^(Q t ,P t ) 

= ^ w (Qt+cft - dtV P G(Q t+dt , P t+dt ), P t+dt + <ftV Q G(Qt+dt, P t+dt )) 
= H^(Q t+dt , P t+dt ) + dt(V P if W • V Q G - V Q tf « • V P G) (7) 



and so, 



an 

= V P F- V Q G- Vq#- V p G = -{H,G} (8) 



which can be re-written as, 

a tt 

- w + {H,G} = (9) 

Note that the evolution equation relates the value of H^ +dt ^> at (Qt+dt, Pt+dt) to the value of 
at (Qt+dt, Pt+dt). Thus, this evolution generates a one-parameter family {H(Q,P;t)} 
of Hamiltonians, whose evolution under the action of the infinitesimal generating functions 
G(Q, P; t) is given by the previous equation (Eq. (9)). 

For what follows in this paper it will prove convenient to represent the dynamics in 
Fourier space. To this end, assume that H is periodic in each Qi with period Lj. Then we 
shall choose G to also be periodic in each Q { with period L { . Define V = L x ■ ■ ■ L D , and let 
Q denote an arbitrary D-dimensional box of side lengths Li, . . . , Lp. Then, 

H(Q,p) = W H ^ p y m]i ' Q ( io ) 



where 



HJP) = / dQH(Q, P) e - 2mk Q (11) 



and similarly for G(Q,P). 

The full evolution done component-wise gives, 



E9Hk 27rik-Q 1 /V"^ ^7 U ^nik'-Q • -ill n 27rik"-Q 

-gf e = 771 (E Vpi/ k /e ^2fi^kG k «e 



1 v ^ 

k at k' k" 

2m k'^e 2mk '- Q • ]T VpG k »e 2mk "- Q ) (12) 



which gives, 

E ^ 2mk Q = ^ E [( k " • VptfkOGk- - (k' • V P G k »)^]e 2 -( k ' +k ")-Q 

97T? 

= -TT EE t( k ' • Vptfk-kOGv " ((k - k') • VpG k 0^k-k']e 27rlk ' Q (13) 

^ k k' 

Our component-wise evolution is therefore, 

—± = 2ni- ]T t(k' • V P // k _ k 0Gk' - ((k - k') • V P G k ,)H k _ k ,\ (14) 

Note that the evolution preserves the integration limits of our original system. Thus the 
topology of the original phase space is preserved. 

As a final note for this subsection, it should be mentioned that degrees of freedom for 
which Li is finite can be treated in an action- angle formalism in which Lj = 1, while degrees 
of freedom for which Lj = oo have their corresponding Fourier sums replaced by integrals. 

2.2 Evolution of the Generating Function 

The above treatment tells us how the Hamiltonian evolves under the action of the infinites- 
imal generating functions, but it tells us nothing of the evolution of S(q, P;t), the overall 
generating function from the initial system (q, p) to (Q t , P t ). We shall deal with this now. 

At time t, our generating function is S(q, P; t), taking us from (q, p) to (Q t , P<). At time 
t + dt, our generating function is S^q, P; t + dt), taking us from (q, p) to (Q t+dt , P t+dt ). We 
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know that (Q t , P t ) and (Qt+dt, Pt+dt) are connected by an infinitesimal generating function 
Q* ■ Pt+dt + dtG(Q t , P t +dt)- We must have, 

P = ^ (q,P ' ;t) = Q^(^ p t + df,t + dt) (15) 

Now, 

dS dS dG 

^ (q,Pt;t) = Q^(^ P t + dt + dt—(Q t ,P t+dt );t) 

dS d 2 S dG 

= dq~ (q ' Pt+dt] ^ + dP9q~ (q ' Pt+dt] ^ ' dQ (Qt ' p ' +d ' )rft 

= dq~ (q ' P * +d * ; ^ + dQ (Q *' P * +d ^ ' dqdP (q ' P * +d * ; t)c?t 

(9S 1 <9G dQ 

= — (q, P t+dt ; + — (Q(q, P t ; t), P t+dt ) • ^(q, P t+dt ; t)dt 

dS dG dO 

= — (q, P t+dt ; t) + — (Q(q, P t+cit ; t), P t+dt ) • -^(q, P t+d t] t)dt 

dS dG 
= —(q,P t+dt ;t) + —(Q(q,P t+dt ;t),P t+dt )dt (16) 

The transformation between the second and third lines in the above derivation was obtained 
by switching from column vectors to row vectors. 
We also have, 

dS dS d 2 S 

— (q, Pt+dt] t + dt) = — (q, P t+dt ; t) + 7^(q, Pt+dt] t)dt (17) 

and so, 

d 2 S dG 



dtdq <9q 

We also have, 



— (Q(q,P t+<ft ;*), Pt+*) (18) 



(9S* dG 

— (q, Pi +d i + eft— (Q t , Pt+dt); t) 

dS d 2 S dG 

— (q, P t+dt ; t) + (q, P t+(it ; t) • — (Q t , Pt+^df (19) 



and 



dS 

Qt+dt = gp(ci:'P t+dt ;t + dt) 



dS d 2 S 

— (q, P t+dt ; t) + (q, Pt+dt; t)dt (20) 



and 



dG 

Qt+dt = Qt+gp(Qt,Pt+dt)dt 

dS d 2 S dG dG 

= gp(<l> p t+df, t) + (q, P t+d u t) ■ gqiQt, Pt+dt)dt + — (Q t , P t+dt )dt (21) 

and so, 

d 2 S d 2 S dG dG 

dtdP t+dt '' ^ ~~ dP 2 t+dt '' ^ ' dQ Pt + dt > + gp^ tj Pt + dt ' 

dG dQ dG 

= — (Q(q, P t ; f), P t+dt ) • ^(q, P t+dt ] t) + ^(Q(q, p t5 *), P *+^) 

<9G <9Q dG 

= ^g(Q(q, P*+dt; £), Pt+dt) • Qpli' Pt+<it! + 0p(Q(q> Pt+dt; t), Pt+dt) 

dG 

= — (Q(q, P t+dt ; t), Pt+dt) I q (22) 

Once again, a switch from column to row vectors gives us the equality between the first and 
second lines of the above derivation. Therefore, we see that, 

d 2 S dG ,dS . „ , 

= ^-(™(q,P;*),P) (23) 



and 



dqdt <9q dP 



d 2 S dG ,dS 



»5(»5(q.P;*),P)l (24) 



so we obtain, 



<9P<% <9P V <9P V ^' ' y ' ' U 



^(q ) P ^) = ^(^(q,P;t),P) + C(t) (25) 
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Since we only care about (Q t , P t ) and not about S itself, we can drop the C(t) term, giving 
us, finally, 

At time t — 0, S is simply the identity transformation, so that S^q, P; 0) = q-P. Thus, S 
is not periodic in each with period Lj. Define nL = (raiLi, . . . , udLd), and note, however, 
that <S(q+ nL, P; 0) = P • nL + S'(q + nL, P; 0). We claim that this property is preserved by 
the evolution. To see this, suppose that at time t we have S^q+nL, P; t) = PnL+S^q, P; t). 
Then, 

— (q + nL,P;t) = nL + — (q,P;t) (27) 
The periodicity of G then implies, 

9S, T „ , OS, n . 

— (q + nL,P;t) = — (q,P;t) (28) 

Therefore, since S^q + nL, P; 0) = P • nL + ^(q, P; 0), it follows from the previous equation 
that this property holds at all t. Thus, although S is not periodic in the we still need 
only track S for q in a D-dimensional box of side lengths Li, . . . , Ld- 



3 Choosing G 

We want an approach that minimizes the dependence of H on Q. The condition that H be 
independent of Q is equivalent to the vanishing of the nonzero Fourier components if k (P). 
We therefore seek to minimize the dependence of H on Q by choosing G in such a way that 
the |if k (P; t)\ are continuously decreasing for all k ^ 0. 
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At some time t, we can write if(Q,P;t) = # (0) (P;t) + # (1) (Q,P;t). The idea is that 
contains the piece of the Hamiltonian which is only dependent on P, and contains 
the remainder. Then, 

= Vp# (0) • V Q G + V P # (1) ■ VqG - V Q # (1) • V P G 

= V P # (0) • V Q G-{H^,G} (29) 

In the limit of a first-order perturbation on an integrable Hamiltonian, the relevant 
equation is, 

™ = V P tf(°) • V Q G (30) 

In Fourier space, this becomes, 

HE = I V ^V**-Q = V P H(°) • ^ V kG k e 2mkQ (31) 



so 



= 27ri(k • Vpif (0) )G k (32) 



9* 

Then = H k ?§± + ^k^r = 27 ™( k ' V F H^)(H k G k - ifkG*). Therefore, in the first- 

order limit, the gradient- descent prescription for minimizing \H k \ 2 = H k H k , k ^ is to set 
G k = 2ni(k ■ V P # (0) )i/ k . Then G k = -2ni(k ■ V P H^)H k , so in the first-order limit we 
obtain that, 

= -8vr 2 (k ■ V F H^ fH k H k (33) 

which is clearly negative. For stronger perturbations, this is no longer the gradient-descent 
prescription. However, for nearly integrable systems (the ones of interest to us in this paper) 
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the perturbation should still be sufficiently weak that the above choice for G k will shrink the 
HkHk, k 7^ 0. Therefore, we take, 

G k = 2m(k ■ V P # (0) )i/ k (34) 

This gives, 

G(Q, P) = ^ E G k (P)e 2 ^ k - Q = ^mV-pH^ • £ kif k e 2 ™ kQ 

= ^V P ff(°» • V Q £ H^ ikQ = Vpff W . Vq if (35) 

so G(Q,P) = V P if (0) • V Q H = V P if (0) • V Q ff (1) . Note that the Fourier expansion of G 
involves terms of the form k- VpH^°\ A k for which k- VpH^ = is a generalized resonance 
at P. The integrability of is destroyed by the resonant terms in H^. 

Our evolution does not formally distinguish between resonances, near-resonances, and 
non-resonances. The evolution is done on the entire Hamiltonian without any terms ne- 
glected. However, the closer a term is to being resonant, the smaller the corresponding 
Fourier component of G, and so the less that term is affected by the evolution. 

In the first-order limit, H ( -°\P; t) differs from if (°)(P; 0) by a correction which is at most 
first-order in Therefore, if we use if (°)(P; 0) instead of if (°)(P; t) in our prescription for 

choosing G in Eq. (35), we get a discrepancy of at most second-order in so that the two 
formulations are equivalent to first-order. Since our prescription for choosing G was derived 
from the first-order limit of the evolution of if, we see that it is equivalent to use if (°)(P; 0) or 
if (°)(P; t) in Eq. (35). Finally, our t — Hamiltonian is usually given as if (q, p) = H (p) + 
y(q>p)> where i?o is the zeroth-order, integrable Hamiltonian, and V is the perturbation. 
We can extract the k = Fourier component of V, writing V(q, p) = V (p) + V"(q, p). 
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Then #(Q,P;0) = H (P) + V (P) + V(Q,P), so that if(°)(P;0) = H (P) + V (P), and 
H^(Q,P) = V"(Q, P). Therefore, note that in the first-order limit it is equivalent to use 
if(°)(P;0) or H (P) in the prescription for choosing G. Once again, this means that it is 
equivalent to use if (°)(P;0) or H (P). In what follows H^(P;t), if(°)(P;0), and if (P) 
will all be denoted by H^°\ or if(°)(P). When required, we will specify to which we are 
referring. 

We conclude this section by deriving the PDE governing the evolution of S, given our 
prescription for choosing G. We have, 
(99 

— (q, P; t) = G(Q(q, P; t), P; t) = V P if <°> • V Q if(Q(q, P; t), P; t) (36) 
Now, we know that H(Q, P; t) = if(q, p; 0) = H(q, §^(q, P; t); 0). We also have, 



d d dQ d d 2 S 



dq dQ dq dQ dqdP 

Switching from row to column vectors gives ^ = dp|q flj' and so, 



(37) 



dQ y dPdq } dq 1 ; 



and so the dynamics of S is governed by, 

dS _ dHW d 2 S dH(q,^;0) 



(39) 



dt dP y dPdq J dq "P 

0H(q "0) II I 2 

Now, l 9q aq ' | p = If | p (q, f§ ; 0) + |f | q (q, §§ ; 0) • f^f. Once again, switching from row 
to column vectors gives, 

d#(q,f;0), dH , , dS n . <9 2 ,S<9fi| , <9S n , 

Note that since if(q, p; 0) is simply our initial Hamiltonian, this PDE for S involves S only. 
We do not need the evolution of H in order to get the evolution of S. 
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The numerical evolution of S is described in Appendix A. The numerical evolution of 
H in the case of action- angle variables is discussed at length in Ref. 1. Since the case for 
arbitrary canonical pairs is handled similarly, we do not give numerical details for the H 
evolution in this paper. 



4 The First-Order Limit 

From our choice of G in the previous section, it follows that in the limit of a first-order 
perturbation on an integrable Hamiltonian, 



= -47r 2 (k- V P # (0) ) 2 # k (41) 



dt 

Our first-order solution yields H k (P; t) = H k (P; 0) exp[-47r 2 (k • V P # (0) ) 2 t]. Note that the 
more non-resonant a term, the faster the exponential decay. In particular, resonances are 
not affected at all. 

We now turn to the evolution of S in the first-order limit. To this end, write S(q, P; t) = 
q • P + G(q, P; t). In what follows we shall work to first-order in G and Then, 

dH I . dS n , dH. _ A . d 2 H . _ n .dG 

flq I p (q ' fl? 0) = 5q" (q ' P; 0) + 5P^ (q ' P; 0) *i (42) 

and 

d 2 S dH ■ ^_. Q) _^_^E\ p dG . Q) _ d 2 G8H 

<9q 2 dp I q q ' 9q ' dq 2 dp I q ^' <9q ' <9q 2 <9P ^' ' 

Now, tf(q,p;0) = iJ( )(p; 0) + H^(q, p; 0). Then §§(q,P;0) = ^(q,P ; 0). To first- 
°rder, **(q,P;0)g[ = S£(q,P;0)g = 0. Finally, to first-order, 0§£ | q (q,P;0) = 
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^■ 9 QP > (P; 0)- Since each of these terms are either first-order in H^> or G, in the first- 
order limit we take ( dpd q ) 1 = 1- Putting everything together gives us that our first-order 
equation is, 

dG _dH® dHW d 2 GdH^ 

~dt~ OP Oq + <9q 2 OP ' ^ ' 

In Fourier space, this becomes, 

^ = 2vu(k • VpH^H^ - 4vr 2 (k • V P tf (0) ) 2 G k (45) 
at 

Since G k (P; 0) = V k, we obtain, 

■tr(l) 

^ ' J 27T(k-Vp//(°)) 1 J 1 j 

Note that G k (P; t) = for all resonant terms. This can be seen by looking at the origi- 
nal ODE from which the solution is derived, or equivalently by noting that lini k . Vp #(o)^ 
(1 - exp[-47r 2 (k • V P # (0) ) 2 t])/(k • V P # (0) ) = 0. We can write, 

This series is convergent, because the exponential term prevents resonances and near reso- 
nances in the denominator from causing the series to diverge. In practice, in the first-order 
limit we only run t out to some finite value, so that the Fourier components of G remain 
small. This smallness criterion depends on what is the maximum allowable G before a first- 
order approximation is no longer deemed accurate, and is therefore dependent on the specific 
error cutoffs for the system at hand. The more non-resonant terms in the Hamiltonian will 
essentially get killed off, while the less non-resonant terms will get reduced somewhat, though 
not completely. Exact resonances will be unaffected by the evolution. In any event we can 
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let t — > oo to get the first-order perturbation theory result, 



G(q,P) 



1 i 
V2n 



E 




27rjk-q 



(48) 



k • Vpif (°) 



where the sum is over all non-resonant k. Note, however, that the convergence of the various 
Fourier components to their t — > oo limits is not uniform, because the time constant for the 
exponential term is proportional to l/(k • Vp-ff"^) 2 . This goes to infinity as k approaches a 
resonance. 

In reality, the above equation must be solved for finite t, and then take the t — > oo limit. 
This prevents any kind of ambiguities in G, something which will be illustrated with one of 
the analytical examples in the next section. 

Before concluding this section, we should note the similarity between this generalized first- 
order perturbation theory and classical first-order perturbation theory in the special case of 
action-angle variables. Indeed, our generalized first-order perturbation theory reduces to 
classical first-order perturbation theory in the case of action-angle variables. One wonders if 
there is a corresponding KAM-like Theorem for more general "tori" than just those described 
by action-angle variables. 

5 Analytical Examples 
5.1 1-D Example 

Consider the case of free propagation perturbed by a weak, localized potential. Our Hamil- 
tonian is, 



H(q,p) 




+ V(q) 



(49) 
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where V(q) is finite in R and decays to as q — > ±00. For above-barrier energies we can 
construct an action function defined everywhere in coordinate space, given by, 

ff( ,,P) = p/>-*!*M)W (50) 

where P = p(±oo). Along a classical trajectory, we have ^ + V(g) = H(q,p) = |^ = so 
the representation of in the corresponding (Q,P) system is simply H(q,p) = H^(P) = 

p2 

2m' 

Because S^g, P) is only determined up to some (possibly P-dependent) constant term, all 
we can uniquely specify is p = || = P(l — 2m ^fa) ) \ . i n the first-order limit, we know from 
the previous section that our PDE method converges to a steady-state. Our method should 
still come close to a steady-state beyond the first-order, yet still weakly perturbed, regime. 
For the strongly perturbed case, we don't know what our method will do. Nevertheless, we 
can show that for all above-barrier energies, our method has a unique steady-state solution, 
and is given by the formula for p(q, P) given above. To prove this, we note that in one 
dimension the evolution equation for S becomes, 

OS P 1 dH . 



dt m dq ' p 

dPdq 1 



(51) 



Setting |j| = 0, and remembering that for energies above the barrier we have P 7^ 0, we get 



the steady-state equation, 

I = (52) 



p 



dq 

Note that this assumes that ^pj^ 7^ 0, 00, something that we will only be able to check once 
we solve our equation. The steady-state equation can be integrated at constant P to give 
H(q, §^) = C(P). Letting q — > ±00 gives us that H becomes |^ = so C{P) = 
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Therefore, + V(q) = £ => f = P(l - 

Note that gpj^ 7^ 0, 00 for all energies above the barrier, so our PDE approach does 
yield the appropriate steady-state solution. Below the barrier, we cannot construct a real 
action function defined for all q, since the momentum p becomes imaginary in the classically 
forbidden region. At a turning point, = 00, so in any event our method for finding the 
steady-state breaks down below the barrier. Thus, both approaches to finding a real action 
function S presented in this section have an equal range of validity, namely, for all energies 
above the barrier. Semiclassically, this means that our PDE approach recovers WKB theory 
for above-barrier energies. Maitra and Heller ]10[ developed a method to compute above- 



barrier reflection coefficients using the WKB wavefunctions as a distorted-wave basis. The 
above derivation shows that this approach is contained as a subcase within our PDE-based 
approach. 

5.2 2-D Example 

Consider the Hamiltonian, 

2 2 

H(x, y, p x , p y ) = + Ae"^ 2 ) (53) 

We wish to construct an action function S(x,y, P x , P y ) in the first-order limit. By the 
symmetry of this problem, we need only consider P y = 0. So, let's write S(x,y, P x , P y ) = 
xP x + yP y + G(x,y, P x , P y ), where G is our first-order correction. We substitute into the 
time-independent Hamilton- Jacobi equation (HJE) to get, 
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J_,p2 + p 2s + + ^y9G + Xe - a {x>W) (54 ) 

2m x y m dx m dy 



Set P y = 0, E = g to get, 



<9G Am 



e -a(* 2 +» 2 ) (55) 



&r P^ 

which can be integrated to give G(x, y, P-r, 0) = — ^e~ ay2 Jf ^ e~ ax ' 2 dx' + C(y, P x ). Now, C 
is determined by the boundary conditions which we impose on this system. The first set of 
boundary conditions we will consider is the requirement that p y = at x = — oo. We have, 



d_G = 2Xma yc _ av , r x + dC 

dy P x J-oo dy 



r x , 2 dC 

/ e~ ax dx' + ^ (56) 
J-oo dy 



Then at x — — oo we get p y = ^ = 0, so C = C(P X ). Since the HJE only depends on 
the momentum field generated by S and not on S itself, a P^-dependent constant term is 
unimportant, so we can set it to 0. Therefore, one first-order solution is, 

Am 



G 1 (x,y,P x ,0) = -^V^f e'^ dx' 

r x J-oo 

Am 2 rV^ x 2 
= __^L e -^ 2 / e -u 2 du ( 57 ) 

P x \/OL J-oo 

An alternative solution is obtained from the boundary condition that p y = —p v 

I x=— oo I x=oo 

Plugging into our expression for ^ gives ^ = — ■pr^ L J^ye~ ay — which can be solved 



for |2 anc i integrated to yield, 



Cfo, P*) = ^|e-^ + 6(P X ) (58) 



As before, we drop the P^-dependent constant term. The result is another first-order solution, 
<Ux,v,P.,0) = -^-^(f^du-f) 

r x yOL J-oo I 

Am 2 / > v / " :E 2 f 2 

= __^ e -^ 2 (/ e~ u du- e~ u du) 

P x \JOL J —oo J —oo 

20 



Figures la,b, plot the momentum fields generated by Gi and G 2 . Note that G\ produces 
the more physically intuitive manifold of trajectories, since they first head straight toward 
the Gaussian bump, and are only deflected as they approach it. The G2 trajectories, in 
contrast, are symmetrical to the left and right of the bump. This is accomplished by initially 
angling the trajectories inward toward the bump. The trajectories curve in and then curve 
away as they reach the bump. 

We now turn to a first-order treatment using the approach derived from our PDE method. 
We begin by using the position formulation of our PDE in the first-order limit, given by Eq. 
(44). We know that the evolution of G goes to a steady-state, so setting ^ = 0, and plugging 
in the terms for our specific potential, yields, 

m [ x + yV) + m 2 1 dx 2 x + dxdy x y + dy 2 y) ~ 1 J 

At P y = we get = ^i xe -^ 2 +y 2 ) . Integrating gives, 

| = -^ + %«) (61) 
If we require p x = P x , then ^ = C(y,P x ) = 0. Our resulting differential 

I x=— 00 ax I x=— 00 

equation for G is therefore identical to the one obtained via the HJE in the first-order limit. 
We now use the Fourier formulation of our first-order, PDE-based approach, by applying 

2 2 ^22 

Eq. (47). The Fourier components of e~ a ^ x +y -* are given by \^e~~^ kx+k y\ Therefore, we 
obtain, 

i\ 7r f°° r°° exp[— — (kl + k 2 )] , .2 (k x p x +k y p y ) 2 

«■ p - ">■• «) = T,a L L x y U\W ^ ~ e -^V-<—» 

x rn tl m 

(62) 
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We now set P y = 0, giving 

G( X ,y,P x At) = r dk x dk y eM -^] 6XP[ ~ - ^ ] (1 - e^W/rn^^y 

ZOti x J — oo J — oo K x 

i\m [a _ av 2 f°° J; expf-^/c 2 ]^ ._^^ t/ m»^fc,x 



2aP x V 7r J-oo k x 

iXm -W r dkx e-^(l - e- 4 ^ p ^ m2 )(2ni f e 2M dx' + -) 
yQ.ll J -co JO k x 



Am pK 



X 



[l e -*y 2 [ x dx > r dkxe -4%^- e -^*Zr£t/™ 2 y*i***' + 
V a jo J-oo 



iXm _ n „2 r°° dk 
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The second integral vanishes, because the integrand is an odd function of k x . This gives us, 

G(x,y,P x ,0;t) = -^L/V"* 2 f tfa' f° dk x e~^(\ - e-^P^yM 

P x V a Jo J-oo 

Am /7T _™,2 r X , , r°° „ , , 7T 2 , 0l r 7T 2 .oiN^fcx' 



V a 



e~ aj/ / dx' I dk x (exp[ /c 2 ]-exp[ 1 k x ])e~ 



Am /^ e - aj/2 
P x \ a 



/ dx'(J— exp[— ax' 2 ] — 

JO V 7T 



x /2 



/ 7r (I + 4P^ /m2) eXp[ i + 4P 2 t/m 2]) (64) 
We can break this integral into two pieces. First let's note the following change of variable: 
Setting u = v^A gives, /* d\^%e~ a}? = ^ J ^ x due~ u2 ' . Therefore, 



G(x, y, P x , 0; t) = --^- e ~ ay2 ( e'^du - f V^ p ^ 2 e^ 2 du) 

P x \/OL Jo Jo 

Am 9 f^ x 



-.e- ay 



r _ / ^ e" u (65) 

P x y^ I m V ' 

v /l+4p2 a t/ m 2 



Letting t — > oo gives us, 



G(x, y, P x , 0) = -—— e -°* \ e~ u du (66) 

Thus, the Fourier space approach yields the G 2 function obtained previously. Note that it 
was necessary to take the t — > oo limit only at the end, in order to avoid any ambiguities in 
evaluating the integrals. 
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6 A Numerical Example 

We chose to test the S evolution numerically on the two-dimensional Pullen- Edmonds Hamil- 
tonian [JJj], given by, 

p 2 + p 2 1 

H(x,y,p x ,p y ) = x y + -(u 2 x x 2 + u 2 y 2 ) + ex 2 y 2 (67) 

This was the two-dimensional system studied numerically in Ref. 1 in testing the GDA 
method. The results there were compared with those obtained in Ref. 2 using the CHM 
method. However, while in Ref. 1 the evolution was done on H in the context of obtaining 
a semiclassical quantum spectrum, in this case it is the PDE for S that is being tested. 

We set u x = uj y = 1. The t = canonical representation is simply the harmonic oscillator 
action-angle basis, denoted by (9 x ,9 y , J x , J y ). The arbitrary t canonical representation is 



denoted by (</) x ,(j) y ,I x ,I y ), so that S = S(9 x ,9 y , I x , I y ;t), with J x = §f, and J, 



as 

y d9 y ■ 



The transformation to the harmonic representation is obtained by setting x = J— cos2tt9 x , 



p x = —^J^sm27i9 x , and similarly for y,p y . The result is, 

H(9 x ,9 y , J x , J y ) = ^-(J !r +J 2/ )+^^(l+cos47r^+cos47r^+^(cos47r(^+^)+cos47r(^-^))) 

(68) 

The only nonvanishing Fourier components are ±(2, 0), ±(0, 2), ±(2, 2), and ±(2, —2). Note 

in particular that ±(2, —2) is a resonance. 

Of the three prescriptions for choosing G, the simplest one to use is — H — 

j-(l, 1) • (I x , I y ), so that Vp-f^ -* is taken to be j-(l, 1). In addition, for this Hamiltonian it 

is readily verified that, 

dH i / dS e dS dS , . t . 1 . . t ,„ . . . „ ,„ . . . 

— q, — ;0 = (sm4 7 rg :c + - sm47r (9 x + 9 y +sm47r (9 x -9 y , 

aq 1 p aq n o9 x o9 y 2 
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sm4ir9y + -(sin 47^ + 9 y ) - sm4n(9 x - 6 V ))) (69) 

and, 

dH 1 e 1 dS dS 

— | q = —(1, 1) + — (l+cos47r^.+cos47r^+-(cos47r(^,+^)+cos47r(^,-^)))(— , — ) 

(70) 

We substitute into Eq. (40) and then into Eq. (39) to get the PDE for S. 

We set e = 0.05, and I x — I y — 9tt. This corresponds to a torus with zeroth-order 
energy E — 9 at t — 0. The system is nearly integrable in this regime @, yet a study 
of the semiclassically obtained energy spectrum [|l], 0] shows clear differences with first- 
order perturbation theory. We therefore test our PDE beyond the first-order limit with this 
example. 

We set N = 20, DX = 0.1, DT = 0.034, and GDSZ = 9, which allowed us to propagate 
the PDE out to a time of 0.306 (the meaning of these parameters is given in Appendix A). 
The rate of change of S on the grid reaches a minimum at this point (which was determined 
by tracking J ((dS/dt) 2 ) on the grid), so the evolution was stopped here. Because our 
PDE approach is only the gradient- descent prescription in the first-order limit, for finite 
perturbations there is no reason to expect the evolution to reach steady-state. We discuss 
this issue further in Appendix C. 

Figures 2a-d show the results of the evolution at times t = 0.0, 0.1, 0.2, and 0.3. It should 
be noted that we are not plotting S in these graphs. Rather, we plot H(6 X) 9 yi Jj^, ■§§-)■ 
There are two reasons for this: First, because we are working in the weakly perturbed, 
though still beyond first-order, regime, S remains fairly close to the identity transformation. 
The effect on H(9 X , 9 y , Jjp, ■§§-), however, is dramatic, and so it is much more convenient to 
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represent the evolution of S in this indirect fashion. Second, even if the perturbation were 
sufficiently strong to significantly deform S, the only way to determine if the deformation of 
S is correct is to look at its effect on H. 

Using the formula derived in Appendix B, it is readily shown that the first-order solution 
to H(9 X , 6 y , gj^, Jj^) is, 

H{e ^W x 'Wj = ^(^ + 4) + ^(l + e- 4t (cos4vr^ + cos4vr^) + 

i(cos47r(^ - 9 V ) + e~ 16t cos4tt(6 x + 9 y ))) (71) 

While the perturbation is sufficiently strong that there are clear differences with first-order 
perturbation theory, the perturbation is still sufficiently weak that the first-order result 
provides a qualitative and semiquantitative picture of how the evolution should proceed. 
Letting t — » oo, we see that the long-time limit of the evolution gives, 

H' I s<'- + + t^ 1 + \ cos ^ - ^ (72) 

This is plotted for our specific set of parameters in Figure 3. Note that the numerical 
evolution does indeed deform S in such a way that the graphs in Figures 2a-d evolve to look 
like the graph in Figure 3. Without the (2, —2) resonance, H is integrable at this energy ||, 
so it is possible to transform to a basis in which H only depends on (I x ,I y ). Instead, the 
evolution eliminates as much of the nonresonant behavior as possible, but the resonant angle 
dependence arising from the (2, —2) term remains. In Figures 4a-d we present the results 
of the evolution on the Hamiltonian obtained by removing the (2, —2) resonance from the 
Pullen-Edmonds term. We changed DT to 0.033 and GDSZ to 11. All other parameters 
are otherwise unchanged. In this case, the evolution should flatten H, and as Figures 4a-d 

25 



confirm, this is exactly what happens. 



7 Conclusions and Future Research 

This paper presented a PDE-based, phase-space deformation approach to optimize the canon- 
ical basis with which to globally represent a nearly integrable Hamiltonian. This paper is a 
generalization of the GDA method in that we are no longer restricted to action-angle vari- 
ables. Any canonical representation of phase space can be used, making the GDA approach 
adaptable to Hamiltonians other than just those describable in an action-angle formalism. 
This paper also extends the GDA method because it now provides an expression for the 
evolution of S, the overall generating function connecting the initial canonical basis to the 
final canonical basis. Because semiclassical wavefunctions are constructed from classical ac- 
tion functions via Eq. (2), we now have an approach for constructing the wavefunctions 
associated with the optimized torus basis. 

The purpose of this paper was to develop the PDE-based approach in complete generality, 
and then to demonstrate it via analytical results in the first-order regime and a simple 
numerical example in the nearly integrable regime. The emphasis in this paper was on the 
PDE for the evolution of S, because the PDE for H was tested extensively in Ref. 1 in the 
case of action-angle variables. 

The PDEs for H and S involve quantities that are purely classical. Thus, there was 
no quantum mechanics in this paper. Nevertheless, as was mentioned previously, the moti- 
vation for this work is derived from semiclassical quantum mechanics. Ref. 1 has already 
used the H evolution as a method to optimize semiclassical torus bases for use in quantum 
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calculations. For future work, we seek to use the S evolution as a tool to actually construct 
the semiclassical wavefunctions. We also hope to apply our approach to real systems. One 
area where this PDE-based approach might be useful is in the determination of vibrational 
energies of polyatomic molecules. 

The generality of the PDE-based approach means that other previous methods are con- 
tained as subcases. As mentioned in Section 5.1, this approach contains as a subcase the 
use of one-dimensional WKB wavefunctions as a distorted-wave basis for computing above- 
barrier reflection coefficients. This method was developed by Maitra and Heller in Ref. 10. 
In their paper, they raised the issue of generalizing their technique to higher dimensions, 
and to action-angle systems. Our PDE-based approach is exactly this generalization, since 
it unifies the Maitra and Heller method and the GDA method under one general approach. 
Furthermore, the PDE-based approach also reduces to classical first-order perturbation the- 
ory in the limit of a first-order perturbation on an integrable Hamiltonian. 

To conclude, we should add that Heller has often made the comparison between the 
separatrix region generated by a local potential bump in one dimension to the resonance zone 



structure in a Poincare surface of section of a nearly integrable Hamiltonian g, [10[ . Indeed, 
Heller regards the above-barrier reflection problem as a prototype for the more complicated 
case of dynamical tunneling between invariant tori facilitated by resonance zones. Because 
both types of systems can be treated within the same PDE-based approach, they are, in 
fact, formally equivalent phenomena. 



27 



8 Acknowldegments 



This research was supported by the Department of Chemistry and Chemical Biology at Har- 
vard, the Department of Physics at Harvard, and by an NSF Graduate Research Fellowship. 
The author would like to thank E.J. Heller for providing the inspiration for this work, and for 
helpful discussions leading to its completion. The author would also like to thank Prof. An- 
thony Yezzi, Jr. (of the Georgia Tech Department of Electrical and Computer Engineering) 
for helpful conversations regarding the numerical solution of the PDEs. 



A Numerical Propagation of S 

If H is given analytically, then the derivatives of H can also be determined analytically. 
Therefore, it can be seen from Eqs. (39) and (40) that the numerical propagation of S only 
requires the numerical evaluation of the partial derivatives of S. This is done using centered 
differences. 

The q-grid is given by {(niLi/N, . . . ,nDLo/N)\ni = 0,...,N — 1}, giving N D grid 
points. Note that the g^'s do not need to go all the way to L i: since S^q + nL, P; t) = 
P ■ nL + S(q, P; t) throughout the evolution. We track all P on a grid of canonical momenta 
about some central momentum P , where our grid consists of all canonical momenta Pk = 
P + DXk, with k = (k 1: . . . , k D ) satisfying \ki\ + . . . + \k D \ < GDSZ. Let us denote this 
set by fi(Po, GDSZ). Since our evolution involves a first-derivative in P of S, we cannot 
compute H at the boundary of the P-grid. The result is that we can only propagate on 
f2(P , GDSZ — 1), so that at each iteration the value of GDSZ shrinks by 1. This collapsing 
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boundary method is described in further detail in Ref. 1, since it also arises naturally in the 
numerical implementation of the if evolution. 

Once ||- has been evaluated on all possible grid points, the propagation by some time step 
DT is done using the Explicit Euler method, which means that we set S(q_, P fc ; t + DT) = 
S(q,P k ;t) + DT§(q,P k ;t). 

Finally, suppose we are considering a system with unbound degrees of freedom, that is, 
some of the = 0. Then we track those qi G {q^ ± nA\n = 0, . . . , Ni}. At each time step, 
we can only compute ^ up to n — iVj — 1, so that after each time step we shrink our set of 
qi by decreasing by one. 

B An Additional First-Order Result 

In this section we derive the first-order result for if(q, §|;0). Following the procedure in 
Section 4, we write S = qP+G(q, P; t). We also write if (q, p; 0) = ff (0) (p; 0)+ffW(q, p; 0). 
Using p = || = P + V q G(q, P; t), we get to first-order that, 

if (q, p; 0) = if (°) (P; 0) + V P ff (0) (P; 0) • V q G(q, P; t) + H « (q, P; 0) (Bl) 

For simplicity, we assume that G was chosen using if^(P) = (P;0). As mentioned 
before, in the first-order limit all three prescriptions for choosing G are equivalent. Using 
the other two prescriptions will lead to at most second-order corrections in our final result. 
From Eq. (47) we get that, 

VpJf <°> • V q G(q, P; t) = -i £ H£\l - e -4- 2 (k-v P ^o) )2t)e2mk . q 

= _|_ ! V ^( 1 ) e -^ 2 ( k -Vp^ (0) ) 2 < e 2T*.q ^ B2 ^ 

^ k^O 
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and so we obtain, 




q,P;t);0) = J ffW( 



P;0) + -E< (P;0) 



e -47r 2 (k-V P H(°)) 2 t e 27rik-q 



(B3) 



C Propagation Time 

Recall from Eq. (32) that the first-order expression for the evolution of H in Fourier space 
is, 



This equation was then used to obtain the gradient-descent prescription for choosing G in 
the first-order limit. For weak perturbations, this prescription no longer coincides with the 
gradient- descent approach, but should still shrink the nonzero Fourier components of H. 
This will occur as long as the right side of Eq. (CI) (or Eq. (32)) is sufficiently dominant 
compared to the remaining terms in the full PDE for H . Note then that for resonances 
and near-resonances this condition does not hold. However, for sufficiently non-resonant 
terms this condition does hold. Thus, in general, for a weak perturbation, our PDE-based 
approach starts out by decreasing the more non-resonant terms of H. The Q-dependence 
of H starts decreasing, and so the rate of change of S decreases as well as the evolution 
proceeds. Eventually, the sufficiently non-resonant terms of H are reduced to a point where 
higher-order terms become important, so that our first-order gradient-descent prescription 
for choosing G will no longer work to reduce the Q-dependence of H. The rate of change 
of S then begins to increase after this point, and eventually the PDE becomes numerically 
unstable. By tracking J {(dS/dt) 2 ) on the grid, it is possible to stop the evolution where the 



2m(k- V P # (0) )G k 



(CI) 



dt 
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rate of change of S reaches its minimum, and consequently where the canonical basis has 
been optimized. 

Of course, the weaker the perturbation, the closer a given k must be to a resonance for 
our choice of G to no longer work to reduce the corresponding if k . Futhermore, the weaker 
the perturbation, the longer it is possible to propagate the PDE before \J {(dS/dt) 2 ) reaches 
its minimum, and the closer this minimum will correspond to a steady-state. It would be 
interesting to develop a simple criterion to estimate at what time this minimum occurs, and 
how far away the system is from steady-state at the minimum. 
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